#include <cmath>
#include <iostream>

#include "gtest/gtest.h"

#include "../RandomRespSolver.h"

// double LogLogIntegral(double x0, double x1, double y0, double y1)
// {
//     double alpha = std::log10(y1 / y0) / std::log10(x1 / x0);

//     if (std::abs(alpha + 1.0) < 1e-8) {
//         return x0 * y0 * std::log(x1 / x0);
//     } else {
//         return y0 / ((1.0 + alpha) * std::pow(x0, alpha)) *
//                (std::pow(x1, alpha + 1.0) - std::pow(x0, alpha + 1.0));
//     }
// }
// SG::Algebra::Real_t LogLogIntegral (SG::Algebra::Real_t x0, SG::Algebra::Real_t x1, SG::Algebra::Real_t y0, SG::Algebra::Real_t y1);
double              RefIntegral (double x0, double x1, double y0, double y1, int N = 1000)
{
    double sum = 0.0;
    double dx  = (x1 - x0) / N;
    for (int i = 0; i < N; i++)
    {
        double xL = x0 + i * dx;
        double xR = x0 + (i + 1) * dx;
        // 对数插值计算每个小段上的 y
        double logyL = std::log (y0) + (std::log (y1) - std::log (y0)) * std::log (xL / x0) / std::log (x1 / x0);
        double logyR = std::log (y0) + (std::log (y1) - std::log (y0)) * std::log (xR / x0) / std::log (x1 / x0);
        double yL    = std::exp (logyL);
        double yR    = std::exp (logyR);
        sum += 0.5 * (yL + yR) * (xR - xL);  // Trapz
    }
    return sum;
}

TEST (LogLogIntegralTest, BasicCases)
{
    // 用手工算的结果验证

    // Case 1: y = x (即 y0/x0 = y1/x1)，alpha=1，测试简单斜率
    {
        double x0 = 1.0, x1 = 10.0;
        double y0 = 1.0, y1 = 10.0;  // y 和 x 成正比

        // 对于 y=x，积分 ∫ y dx = ∫ x dx = (x^2)/2, 结果是 (100-1)/2=49.5
        double expected = (std::pow (x1, 2) - std::pow (x0, 2)) / 2.0;
        double result   = RandomRespSolver::LogLogIntegral (x0, x1, y0, y1);
        EXPECT_NEAR (result, expected, 1e-10);
    }

    // Case 2: 常数函数 y=2 (y0=y1)，应该是 2*(x1-x0)
    {
        double x0 = 1.0, x1 = 5.0;
        double y0 = 2.0, y1 = 2.0;
        double expected = 2.0 * (x1 - x0);
        double result   = RandomRespSolver::LogLogIntegral (x0, x1, y0, y1);
        EXPECT_NEAR (result, expected, 1e-10);
    }

    // Case 4: 特殊情况 alpha ≈ -1
    {
        double x0 = 1.0, x1 = 10.0;
        double y0 = 10.0, y1 = 1.0e-1;  // y1/y0 和 x1/x0 设计成 alpha ≈ -1
        double result = RandomRespSolver::LogLogIntegral (x0, x1, y0, y1);
        // 这里只检查结果不为NaN或inf
        EXPECT_TRUE (std::isfinite (result));
    }
}

TEST (LogLogIntegralTest, DecreasingY)
{
    double x0 = 1.0, x1 = 10.0;
    double y0 = 10.0, y1 = 1.0;

    double expected = RefIntegral (x0, x1, y0, y1);   // 高精度参考值
    double result   = RandomRespSolver::LogLogIntegral (x0, x1, y0, y1);
    EXPECT_NEAR (result, expected, expected * 1e-3);  // 误差控制在0.1%
}